FOV Registration over Multiple Days¶
Progress Report Sept. 2019
In [37]:
plt.figure(figsize=[5,5])
# r, g, and b are 512x512 float arrays with values >= 0 and < 1.
from PIL import Image
rChannel = out_previous[-3]
gChanel = out_previous[-2]
bChannel = out_previous[-1]
rgbArray = np.zeros((315,315,3), 'uint8')
rgbArray[..., 0] = rChannel*(rChannel>0.64)*256
rgbArray[..., 1] = gChanel*(gChanel>0.64)*256
rgbArray[..., 2] = bChannel*(bChannel>0.64)*256
img = Image.fromarray(rgbArray)
plt.imshow(rgbArray)
#plt.imshow(,alpha=0.1,cmap='Reds')
Out[37]:
In [464]:
f,(ax,bx) = plt.subplots(1,2,figsize=[40,80],sharex=True,sharey=True)
A = cnmr.get_A_sparse(resultList[0])
contours = cnmr.get_contours(A,nrgthr=0.75)
ax.imshow(stack[0])
#bx.imshow(stack[0].T)
for cont in contours:
ax.plot(cont[:,0],cont[:,1],'k',lw=1)
cm.utils.visualization.plot_contours(A,stack[0],thr_method='nrg',maxthr=0.32,nrgthr=0.5,ax=bx);
#bx.set_xlim(100,150)
In [448]:
f,(ax,bx,cx) = plt.subplots(1,3,figsize=[40,60],sharex=True,sharey=True)
sr = pystackreg.StackReg(pystackreg.StackReg.RIGID_BODY)
# register each frame to the previous (already registered) one
# this is what the original StackReg ImageJ plugin uses
TrMatrix = sr.register_stack(stack[5:7], reference='first')
out_previous = sr.register_transform_stack(stack[5:7], reference='first')
################################################
RefIdx = 0
bkgIdx = 1
ARef = cnmr.get_A_sparse(resultList[5])
contours = cnmr.get_contours(ARef,nrgthr=0.85)
cont = contours[1254]
cellContour = np.append(cont,np.zeros([cont.shape[0],1]),1).reshape(-1,3)
######################################################
ax.imshow(out_previous[RefIdx],cmap='Greys_r',)
bx.imshow(out_previous[bkgIdx],cmap='Greys_r',)
for _ii,Mx in enumerate(TrMatrix[:]):
newCoord = np.dot(cellContour,TrMatrix[0])
for _jj in range(_ii):
newCoord = np.dot(newCoord,TrMatrix[_jj])
newCont = np.dot(cellContour,Mx)
#print(newCoord)
#plt.plot(point[0],point[1],'k.')
#bx.plot(newCoord[:,0],newCoord[:,1],lw=1)
if _ii==bkgIdx:
bx.plot(newCoord[:,0],newCoord[:,1],lw=2,c='r')
bx.plot(newCont[:,0],newCont[:,1],'k',lw=1)
if _ii==RefIdx:
ax.plot(cont[:,0],cont[:,1],lw=2,c='r')
ax.plot(newCoord[:,0],newCoord[:,1],lw=2,c='k')
#ax.set_xlim(100,150)
#ax.set_ylim(250,315)
rChannel = out_previous[RefIdx]
gChanel = out_previous[bkgIdx]
#bChannel = out_previous[-1]
rgbArray = np.zeros((315,315,3), 'uint8')
rgbArray[..., 0] = rChannel*(rChannel>0.4)*256
rgbArray[..., 1] = gChanel*(gChanel>0.4)*256
#rgbArray[..., 2] = bChannel*(bChannel>0.64)*256
img = Image.fromarray(rgbArray)
cx.imshow(rgbArray)
Out[448]:
In [590]:
sr = pystackreg.StackReg(pystackreg.StackReg.RIGID_BODY)
f,axarr = plt.subplots(3,(len(stack)-1),figsize=[80,21],sharex=True,sharey=True)
for _ii,item in enumerate(stack[:-1]):
out_previous = sr.register_transform_stack(stack[_ii:_ii+2], reference='first')
ax = axarr[0,_ii]
bx = axarr[1,_ii]
cx = axarr[2,_ii]
#print(len(out_previous))
ax.imshow(out_previous[0])
bx.imshow(out_previous[1])
rChannel = out_previous[0]
gChanel = out_previous[1]
rgbArray = np.zeros((315,315,3), 'uint8')
rgbArray[..., 0] = rChannel*(rChannel>0.4)*256
rgbArray[..., 1] = gChanel*(gChanel>0.4)*256
img = Image.fromarray(rgbArray)
cx.imshow(rgbArray)
ax.set_aspect(1.0)
#axarr[1,6].set_aspect(1.0)
In [588]:
sr = pystackreg.StackReg(pystackreg.StackReg.RIGID_BODY)
RegStack = sr.register_transform_stack(stack, reference='previous')
f,axarr = plt.subplots(2,(len(stack)-1),figsize=[80,14],sharex=True,sharey=True)
for _ii,item in enumerate(stack[:-1]):
out_previous = sr.register_transform_stack(stack[_ii:_ii+2], reference='first')
ax = axarr[0,_ii]
bx = axarr[1,_ii]
rChannel = out_previous[0]
gChanel = out_previous[1]
rgbArray = np.zeros((315,315,3), 'uint8')
rgbArray[..., 0] = rChannel*(rChannel>0.4)*256
rgbArray[..., 1] = gChanel*(gChanel>0.4)*256
img = Image.fromarray(rgbArray)
ax.imshow(rgbArray)
###############################
rChannel = RegStack[_ii]
gChanel = RegStack[_ii+1]
rgbArray = np.zeros((315,315,3), 'uint8')
rgbArray[..., 0] = rChannel*(rChannel>0.4)*256
rgbArray[..., 1] = gChanel*(gChanel>0.4)*256
img = Image.fromarray(rgbArray)
bx.imshow(rgbArray)
ax.set_aspect(1.0)
#axarr[1,6].set_aspect(1.0)
In [592]:
baseFolder = '/media/alireza_chenani/dataWork/stress_project/315_stress_week/BL6-111/'
tiffList = sorted([item for item in locate('C*.tif',baseFolder)])
stack = np.array([io.imread(tif).T for tif in tiffList])
RegStack = sr.register_transform_stack(stack, reference='previous')
Randsin = np.random.rand(500)/2
for _ii,img in enumerate(RegStack):
prevImg = RegStack[_ii-1]
temp1 = prevImg*(prevImg>np.quantile(prevImg,0.9))
temp2 = img*(img>np.quantile(img,0.9))
overlap = np.nansum((temp1*temp2))/np.nansum((temp1*temp1))#np.nansum((temp1*temp2)/(temp1*temp1))
ov = []
for sin in Randsin:
R = np.array([
[1-sin**2, -1*sin, 0],
[sin, 1-sin**2, 0],
[0, 0, 1]])
Rimg = sr.transform(prevImg,tmat=R)
tempR = Rimg*(Rimg>np.quantile(prevImg,0.9))
ov.append(np.nansum((temp1*tempR))/np.nansum((temp1*temp1)))
Thr = np.quantile(ov,0.95)
if _ii<8:
CLR = 'k'
else:
CLR = 'r'
plt.scatter(_ii,overlap/Thr,c=CLR)
#plt.xlim(.5,14.5)
#plt.ylim(.9,8012.5)
plt.axhline(1)
Out[592]:
In [587]:
Randsin = np.random.rand(500)/2
img = stack[0]
ov= []
for sin in Randsin:
R = np.array([
[1-sin**2, -1*sin, 0],
[sin, 1-sin**2, 0],
[0, 0, 1]])
Rimg = sr.transform(img,tmat=R)
temp2 = Rimg*(Rimg>np.quantile(img,0.9))
temp1 = img*(img>np.quantile(img,0.9))
overlap = np.nansum((temp1*temp2))/np.nansum((temp1*temp1))
ov.append(overlap)
print(np.quantile(ov,0.95))
In [550]:
sns.distplot(ov)
plt.axvline(np.quantile(ov,0.9))
Out[550]:
In [548]:
np.nansum((temp1*temp2))/np.nansum((temp1*temp1))
Out[548]:
In [546]:
315*315
Out[546]:
In [ ]: